#set working directory, clear workspace
  #setwd("")
  rm(list=ls())

#installing/loading necessary pacakges
#install.packages("R2WinBUGS")
#install.packages("foreign")
#install.packages("mcmcplots")
  library(coda)
  library(foreign)
  library(R2WinBUGS)
  library(mcmcplots)
  library(plyr)

##################################################################################################################################  

#reading in matrix of votes and economic variables
  dat <- read.csv("hungary-vote-data-final.csv")
  dat <- dat[,1:12]
  
#defining the inflation gap as difference between two year forecast and current target, deleting observations with NA
  infgap <- dat$twoyearfc-dat$target
  dat <- cbind(dat,infgap)
  dat <- na.omit(dat)

#defining variable names to match BUGS file
  y <- dat$vote
  nvote <- length(y)
  bank <- dat$cbid
  nbank <- max(bank)
  sig <- dat$vote.unc
  lag <- dat$currentrate
  ygap <- dat$ygap
  igap <- dat$infgap
  xr <- dat$xr

#reading in the model file
  model.file <- "rcm-corr.txt"

#setting intial values
  inits <- function() { list(a=runif(nbank,-1,1),
                           bsig=runif(nbank,-1,1),
                           blag=runif(nbank,-1,1),
                           bygap=rnorm(1),
                           bigap=rnorm(1),
                           bxr=rnorm(1),
                           
                           sigma.a=runif(1),
                           sigma.bsig=runif(1),
                           rho=runif(1),
                           sigma.blag=runif(1),
                           sigma.bygap=runif(1),
                           sigma.bigap=runif(1),
                           sigma.bxr=runif(1))  }

#running model in WinBUGS
  bugs.model <- bugs(data=c("y","nvote","bank","nbank", "sig","lag","ygap","igap","xr"),
                    inits=inits,
                    parameters.to.save=c("a","bsig","blag","bygap","bigap","bxr",
                                         "mu.a","mu.bsig","mu.blag","mu.bygap","mu.bigap","mu.bxr",
                                         "sigma.a","sigma.bsig","sigma.blag","sigma.bygap","sigma.bigap","sigma.bxr","rho"),
                    model.file= model.file,
                    n.chains=3,
                    n.iter=150000,
                    n.burnin=50000,
                    n.thin=5,
                    debug=TRUE)

#Extracting Coefficients from bugs.model for Table 1:
  bugs.model #Deviance, DIC
  bugs.model$mean$mu.a
  bugs.model$sd$mu.a
  bugs.model$mean$mu.bsig
  bugs.model$sd$mu.bsig
  bugs.model$mean$mu.blag
  bugs.model$sd$mu.blag
  bugs.model$mean$bygap
  bugs.model$sd$bygap
  bugs.model$mean$bigap
  bugs.model$sd$bigap
  bugs.model$mean$bxr
  bugs.model$sd$bxr
  
#Checking for Credible Intervals on Individual Level "bsig" Estimates
  #For 95% Credible Interval
    
    names <- c("P. Adamecz","H. Auth","A. Balog","T. Banfi","A. Bartfai Mager","P. Bihari","V. Bihari","J. Cinkotai","C. Csaki","F. Gerhardt",
             "I. Hardy","Z. Jarai","B. Kadar","C. Kandracs","F. Karvalits","J. Kiraly","G. Kocziszky","G. Kopits","G. Matolcsy","J. Nemenyi",
             "G. Oblath","G. Pleschinger","A. Simor","G. Szapary","L. Windisch","G. Bager",
             "M. Nagy")
    par(mar=c(3,5.5,1,1))
    caterplot(mcmcout=bugs.model,"bsig",quantiles=list(outer=c(0.01,0.99),inner=c(0.025,0.975)),
              labels=names,col="black",style="plain")
    abline(v=0,lty=2)
  
#Saving Individual Preference Estimates --- Creates File "Hungary_Pref.csv"
  b.int <- bugs.model$mean$a
  b.sig <- bugs.model$mean$bsig
  b.lag <- bugs.model$mean$blag
  b.igap <- rep(bugs.model$mean$bigap,27)
  b.xr <- rep(bugs.model$mean$bxr,27)
  b.ygap <- rep(bugs.model$mean$bygap,27)
  
  names <- c("P. Adamecz","H. Auth","A. Balog","T. Banfi","A. Bartfai Mager","P. Bihari","V. Bihari","J. Cinkotai","C. Csaki","F. Gerhardt",
             "I. Hardy","Z. Jarai","B. Kadar","C. Kandracs","F. Karvalits","J. Kiraly","G. Kocziszky","G. Kopits","G. Matolcsy","J. Nemenyi",
             "G. Oblath","G. Pleschinger","A. Simor","G. Szapary","L. Windisch","G. Bager",
             "M. Nagy")
  
  pref.parameters <- cbind(names,b.int,b.sig,b.lag,b.igap,b.xr,b.ygap)
  
#Generating fitted values for individual central bankers --- Used in Figure #5
  cap.R <- b.int + (b.sig*(mean(dat$vote.unc))) + (b.lag*(mean(dat$currentrate))) + (b.igap*(mean(dat$infgap))) + (b.xr*(mean(dat$xr))) + (b.ygap*(mean(dat$ygap)))

#Adding Information about Appointing Prime Minister and the Level of Uncertainty At the Time of Appointment
  appoint.dat <- read.csv("appointunc.csv")
  pmcode <- appoint.dat$pmcode
  appoint.unc <- appoint.dat$appoint.unc
  appoint.year <- appoint.dat$appoint.year
  pref.parameters <- cbind(pref.parameters,cap.R,pmcode,appoint.unc,appoint.year)                 

#Write CSV file with individual preferences and appointment info
  write.csv(pref.parameters,file="Hungary_Pref.csv")

##################################################################################################################################
  
dat.pref <- read.csv("Hungary_Pref.csv")
  
#Estimate Models from Table 2
  model1 <- lm(dat.pref$b.int~dat.pref$appoint.unc)
  summary(model1)
  
  model2 <- lm(dat.pref$b.int~dat.pref$appoint.unc+factor(dat.pref$pmcode)-1)
  summary(model2)
    #pmcode3 = Orban
    #pmcode2 = Gyurcsany
    #pmcode1 = Medgyessy
  
##################################################################################################################################

#Robustness of results dropping appointments once there were widespread concerns of corruption beginning with 2013 appointment
  #of Matolcsy
  
dat.pref.early <- subset(dat.pref,dat.pref$appoint.year < 2013)
  
model1.early <- lm(dat.pref.early$b.int~dat.pref.early$appoint.unc)
summary(model1.early)

model2.early <- lm(dat.pref.early$b.int~dat.pref.early$appoint.unc+factor(dat.pref.early$pmcode)-1)
summary(model2.early) 
  
##################################################################################################################################

#Distribution of Votes Cast in ``Ineffective Region''

dat <- read.csv("hungary-vote-data-final.csv")
dat <- data.frame(cbind(dat$meetid, dat$vote.unc, dat$appoint.unc))
dat <- na.omit(dat)
  colnames(dat)<- c("meetid","vote.unc","appoint.unc")
  diff <- dat$appoint.unc-dat$vote.unc
  dat <- cbind(dat,diff)

id <- c(seq(1,37),seq(39,131))
  dat <- dat[order(dat$meetid),]
  vote.count <- count(dat,"meetid")

total.votes <- NULL
  for(i in 1:130){
    seq <- rep(x=vote.count$freq[i],length.out=vote.count$freq[i])
    total.votes <- c(total.votes,seq)
  }
  dat <- cbind(dat,total.votes)

inef <- NULL
  for(i in id){
    count <- sum(dat$diff >= .2 & dat$meetid==i)
    seq <- rep(count,length.out=sum(dat$meetid==i))
    inef <- c(inef,seq)
  }
  dat <- cbind(dat,inef)

dat$maj <- dat$inef/dat$total.votes
  maj.meeting <- aggregate(dat$maj,list(dat$meetid),mean)
  sum(maj.meeting[,2] >= .5)/130

sum(dat$diff >=.2)/nrow(dat)


